O objetivo deste trabalho é realizar uma análise minuciosa em dados imobiliários da cidade de São Paulo no ano de 2017. Através da limpeza e a vizualização dos relacionamentos entre as variáveis disponíves, a ideia é tentar captar insights importantes e, em seguida, construir um modelo de precificação baseado nas informações encontradas.
Para iniciar a análise, começarei pela abertura do arquivo através da função base do sistema e, em seguida, com o auxílio do pacote dplyr, tranformá-la-ei em um tbl.df para obter as melhores possibilidades de manipulação de dados.
# Carrehando pacote necessario
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
# Obtencao dos dados
imoveis <- as.tbl(read.csv2('Crawler_vivareal_venda_2017_desafio.csv', stringsAsFactors = F, encoding = "UTF-8",
na.strings = ""))
glimpse(imoveis)
## Observations: 7,296
## Variables: 19
## $ X.U.FEFF.date <chr> "25/07/17", "25/07/17", "25/07/17", "25/07/17"...
## $ rua <chr> "Coronel Artude Paula Ferreira", "Coronel Artu...
## $ numero <chr> NA, "95", "148", "272", "450", "506", "554", "...
## $ bairro <chr> "Vila Nova Conceição", "Vila Nova Conceição", ...
## $ condominio <chr> NA, "Condomínio Loose In Vila Nova ", "Condomí...
## $ price <int> 1000000, 1000000, 1000000, 1000000, 1000000, 1...
## $ area <int> 71, 67, 45, 79, 70, 70, 48, 85, 90, 100, 73, 7...
## $ price_by_sqm <chr> "14.084", "14.925", "22.222", "12.658", "14.28...
## $ condominium_fee <chr> "R$1.450 ", "R$1.280 ", "R$700 ", "R$648 ", "R...
## $ iptu <chr> "R$280 ", "R$2.200 ", NA, "R$220 ", "R$240 ", ...
## $ rooms <int> 2, 2, 1, 1, 1, 1, 1, 1, 1, 3, 2, 1, 2, 2, 3, 1...
## $ bathrooms <int> 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 2, 2, 2, 2, 4, 1...
## $ garage <int> 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1...
## $ agent <chr> "Riala Empreendimentos Imobiliários", "Nalva I...
## $ agent_number <chr> "(11) 5574-0888 ", "(11) 95772-7990 ", "(11) 9...
## $ latitude <chr> "-235.969", "-23.596.244", "-2.359.186", "-2.3...
## $ longitude <chr> "-466.699", "-46.669.762", "-46.679.232", "-46...
## $ property_d <int> 74376409, 78347836, 77289482, 73325473, 728155...
## $ url <chr> "https://www.vivareal.com.br/imovel/apartament...
O primeiro passo será a limpeza dos dados. Após realizar a abertura do arquivo, identifiquei as necessidades de correção nas variáveis, das quais se destacam: problema nos encodings das variáveis data, latitude e longitude; e o formato de texto das variáveis numéricas. Os pacotes selecionados para tais tarefas foram o lubridate e stringr. Ao término do processo, percebe-se a diferença no novo dataframe.
### Data celaning
imoveis_clean = imoveis
# Pacotes necessários
library(stringr)
library(lubridate)
##
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
##
## date
# data
imoveis_clean <- imoveis
# data
names(imoveis_clean) <- str_replace(names(imoveis_clean), 'X.U.FEFF.date', 'Data')
imoveis_clean$Data <- as.Date(imoveis$X.U.FEFF.date, format = '%d/%m/%Y')
year(imoveis_clean$Data) <- year(imoveis_clean$Data)+2000
# lat e long
imoveis_clean$longitude <- gsub(".", "", imoveis$longitude, fixed = TRUE) %>%
sub("(\\d{2})", "\\1\\.", .) %>%
as.numeric()
imoveis_clean$latitude <- gsub(".", "", imoveis$latitude, fixed = TRUE) %>%
sub("(\\d{2})", "\\1\\.", .) %>%
as.numeric()
# price_by_sqm
imoveis_clean$price_by_sqm <- as.numeric(imoveis_clean$price_by_sqm)
# condominium_fee e iptu
imoveis_clean$condominium_fee <- str_replace_all(imoveis$condominium_fee, 'R\\$', '') %>%
str_replace_all('[:punct:]', '') %>%
as.numeric()
imoveis_clean$iptu <- str_replace_all(imoveis$iptu, 'R\\$', '') %>%
str_replace_all('[:punct:]', '') %>%
as.numeric()
# visualizacao do novo df
glimpse(imoveis_clean)
## Observations: 7,296
## Variables: 19
## $ Data <date> 2017-07-25, 2017-07-25, 2017-07-25, 2017-07-2...
## $ rua <chr> "Coronel Artude Paula Ferreira", "Coronel Artu...
## $ numero <chr> NA, "95", "148", "272", "450", "506", "554", "...
## $ bairro <chr> "Vila Nova Conceição", "Vila Nova Conceição", ...
## $ condominio <chr> NA, "Condomínio Loose In Vila Nova ", "Condomí...
## $ price <int> 1000000, 1000000, 1000000, 1000000, 1000000, 1...
## $ area <int> 71, 67, 45, 79, 70, 70, 48, 85, 90, 100, 73, 7...
## $ price_by_sqm <dbl> 14.084, 14.925, 22.222, 12.658, 14.285, 14.285...
## $ condominium_fee <dbl> 1450, 1280, 700, 648, 1090, 1100, 1100, 800, 7...
## $ iptu <dbl> 280, 2200, NA, 220, 240, 240, 120, 140, 100, N...
## $ rooms <int> 2, 2, 1, 1, 1, 1, 1, 1, 1, 3, 2, 1, 2, 2, 3, 1...
## $ bathrooms <int> 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 2, 2, 2, 2, 4, 1...
## $ garage <int> 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1...
## $ agent <chr> "Riala Empreendimentos Imobiliários", "Nalva I...
## $ agent_number <chr> "(11) 5574-0888 ", "(11) 95772-7990 ", "(11) 9...
## $ latitude <dbl> -23.59690, -23.59624, -23.59186, -23.59758, -2...
## $ longitude <dbl> -46.66990, -46.66976, -46.67923, -46.67164, -4...
## $ property_d <int> 74376409, 78347836, 77289482, 73325473, 728155...
## $ url <chr> "https://www.vivareal.com.br/imovel/apartament...
summary(imoveis_clean)
## Data rua numero
## Min. :2017-07-25 Length:7296 Length:7296
## 1st Qu.:2017-07-25 Class :character Class :character
## Median :2017-07-25 Mode :character Mode :character
## Mean :2017-07-25
## 3rd Qu.:2017-07-25
## Max. :2017-07-25
##
## bairro condominio price
## Length:7296 Length:7296 Min. : 18000
## Class :character Class :character 1st Qu.: 1200000
## Mode :character Mode :character Median : 2900000
## Mean : 4253352
## 3rd Qu.: 5624205
## Max. :48000000
##
## area price_by_sqm condominium_fee iptu
## Min. : 18.0 Min. : 2.432 Min. : 1 Min. : 1
## 1st Qu.: 86.0 1st Qu.: 13.571 1st Qu.: 1071 1st Qu.: 262
## Median : 166.0 Median : 17.026 Median : 2000 Median : 800
## Mean : 215.3 Mean : 17.709 Mean : 2860 Mean : 2074
## 3rd Qu.: 293.0 3rd Qu.: 20.755 3rd Qu.: 4000 3rd Qu.: 1887
## Max. :30700.0 Max. :224.000 Max. :460000 Max. :50000
## NA's :40 NA's :40 NA's :934 NA's :1865
## rooms bathrooms garage agent
## Min. : 1.00 Min. : 1.000 Min. : 1.000 Length:7296
## 1st Qu.: 2.00 1st Qu.: 2.000 1st Qu.: 1.250 Class :character
## Median : 3.00 Median : 3.000 Median : 3.000 Mode :character
## Mean : 2.91 Mean : 3.451 Mean : 2.991
## 3rd Qu.: 4.00 3rd Qu.: 5.000 3rd Qu.: 4.000
## Max. :10.00 Max. :54.000 Max. :16.000
## NA's :16 NA's :761 NA's :194
## agent_number latitude longitude property_d
## Length:7296 Min. :-23.64 Min. :-46.70 Min. :26572520
## Class :character 1st Qu.:-23.60 1st Qu.:-46.67 1st Qu.:68891954
## Mode :character Median :-23.59 Median :-46.67 Median :76929697
## Mean :-23.59 Mean :-46.67 Mean :73681658
## 3rd Qu.:-23.59 3rd Qu.:-46.67 3rd Qu.:81281218
## Max. :-23.56 Max. :-46.64 Max. :83874638
## NA's :1661 NA's :1661
## url
## Length:7296
## Class :character
## Mode :character
##
##
##
##
Com a preparacao dos dados concluída, pode-se prosseguir à análise exploratória. Para tal objetivo, procederemos a uma série de visualizações para ilustrar melhor a relação entre a variável de interesse em relação as demais.
## Visualizacao interativa
library(leaflet)
leaflet(data = imoveis_clean) %>% addTiles() %>%
addMarkers(~longitude, ~latitude, popup = ~as.character(property_d),
clusterOptions = markerClusterOptions()
)
## Warning in validateCoords(lng, lat, funcName): Data contains 1661 rows with
## either missing or invalid lat/lon values and will be ignored
A começar com uma visualização interativa dos imóveis no mapa da cidade, através do pacote leaflet, consegue-se obersvar uma série de 20 coordenadas geográficas distantes da maioria das observações, configurando um ponto de interesse para se analisar com mais calma em seguida.
# visualizacao das variaveis continuas
library(ggplot2)
cont_vars <- c('area', 'iptu','condominium_fee')
plots_cont <- list()
for (var in cont_vars) {
plots_cont[[var]] <- ggplot(imoveis_clean, aes_string(log(imoveis_clean[[var]]), log(imoveis_clean[['price']]))) +
geom_point() +
geom_smooth(method = lm, se = F) +
ggtitle(paste(var, 'x price')) +
theme_minimal()
print(plots_cont[[var]])
}
## Warning: Removed 40 rows containing non-finite values (stat_smooth).
## Warning: Removed 40 rows containing missing values (geom_point).
## Warning: Removed 1865 rows containing non-finite values (stat_smooth).
## Warning: Removed 1865 rows containing missing values (geom_point).
## Warning: Removed 934 rows containing non-finite values (stat_smooth).
## Warning: Removed 934 rows containing missing values (geom_point).
Com o auxílio do pacote ggplot2, pode-se vislumbrar a relação gráfica entre as variáveis contínuas do dataframe. É interessante notar a ocorrência de outliers na relação entre ‘price’, ‘area’ e ‘condominium_fee’, gerando outras observações a serem olhadas mais de perto.
# visualizacao das variaveis factor
fact_vars <- c('bairro', 'rooms', 'bathrooms','garage')
plots_fact <- list()
for (var in fact_vars) {
plots_fact[[var]] <- ggplot(imoveis_clean, aes_string(x = as.factor(imoveis_clean[[var]]), log(imoveis_clean[['price']]))) +
geom_boxplot() +
ggtitle(paste(var, 'x price')) +
theme_minimal()
print(plots_fact[[var]])
}
plots_fact_c <- list()
for (var in fact_vars) {
plots_fact_c[[var]] <- ggplot(imoveis_clean, aes_string(x = as.factor(imoveis_clean[[var]]))) +
geom_bar(stat = 'count') +
ggtitle(paste('Total Observations in ',var)) +
theme_minimal()
print(plots_fact_c[[var]])
}
Por fim toma-se nota das relações entre as variáveis categóricas com auxílio do boxplots e gráficos de barra. Praticamente, todas as variáveis possuem algum outlier interessante, ou poucas observações de valores muito elevados, assim se fará uma sensível análise caso a caso para entender melhor suas particularidades. Outro ponto importante é a inexistência de outros bairros a não ser Vila Nova Conceição.
Nesta seção, realizaremos um breve estudo dos outliers identificados ao longo da análise exploratória para melhor entender seus aparecimentos.
## Investigacao outliers
imoveis_clean %>% filter(log(area) > 10 | log(price) < 10 | log(condominium_fee) > 10)
## # A tibble: 3 x 19
## Data rua numero bairro condominio price area price_by_sqm
## <date> <chr> <chr> <chr> <chr> <int> <int> <dbl>
## 1 2017-07-25 <NA> <NA> Vila ~ <NA> 6.90e6 30700 224
## 2 2017-07-25 Coro~ <NA> Vila ~ <NA> 1.80e4 226 79
## 3 2017-07-25 Mont~ <NA> Vila ~ "Condomín~ 4.80e5 32 15
## # ... with 11 more variables: condominium_fee <dbl>, iptu <dbl>,
## # rooms <int>, bathrooms <int>, garage <int>, agent <chr>,
## # agent_number <chr>, latitude <dbl>, longitude <dbl>, property_d <int>,
## # url <chr>
imoveis_clean %>% filter(bathrooms == 54)
## # A tibble: 1 x 19
## Data rua numero bairro condominio price area price_by_sqm
## <date> <chr> <chr> <chr> <chr> <int> <int> <dbl>
## 1 2017-07-25 Balt~ <NA> Vila ~ "Condomín~ 3.18e6 220 14.5
## # ... with 11 more variables: condominium_fee <dbl>, iptu <dbl>,
## # rooms <int>, bathrooms <int>, garage <int>, agent <chr>,
## # agent_number <chr>, latitude <dbl>, longitude <dbl>, property_d <int>,
## # url <chr>
# outliers de area e price: informações muito discrepantes em relação à realidade. Provavelmente, erros de digitação
imoveis_clean$area[which(imoveis$area == 30700)] <- NA
imoveis_clean <- imoveis_clean %>% filter(!(price == 18000))
imoveis_clean$condominium_fee[which(imoveis$condominium_fee == 460000)] <- NA
# o apartamento com mais de 54 banheiros parece muito fora da media em relacao ao seu preco e area
imoveis_clean$bathrooms[which(imoveis$bathrooms == 54)] <- NA
Os valores de ‘price’, ‘area’ e ‘condominium_fee’, assim como vizualizados acima, estavam muito discrepantes da realidade sem motivo aparente, dando uma boa evidência de serem erros de digitação. Por não haver forma melhor de checar sua veracidade, e para não gerar viés indesejado na análise, optei por retirá-los da amostra. Como o único imóvel com um número muito acima da média de banheiros não seguia outros em relação a área e valor, optei pela mesma abordagem também.
Para a construção de um modelo de precificação eficiente, tem-se de extrair o máximo de informação de cada uma das variáveis e, ao mesmo tempo, evitar o overfitting, portanto a construção de novas variáveis é um passo fundamental da análise.
A começar por features já existentes, faz-se a agregação de algumas variáveis categóricas com o intuito de diminuir classes com pouca representatividade na amostra como visto nos gráficos de barra. A ideia é evitar o overfitting.
# agregação dos valores para não causar overfitting pela baixa quantidade de valores
roomsBathGarage <- imoveis_clean %>%
select(c(rooms, bathrooms, garage, property_d)) %>%
mutate(rooms_coerc = ifelse(rooms > 4, 4, rooms),
bathrooms_coerc = ifelse(bathrooms > 6, 6, bathrooms),
garage_coerc = ifelse(garage > 5, 5, garage)) %>%
select(-c(rooms, bathrooms, garage))
Como as variáveis de identificação, em seu estado bruto, não são relevantes à análise, executarei a construção de outras com possibilidade maior de fornecer informação. Após realizar sua criação, faço plots demonstrando suas informações em relação às variáveis antigas.
# criacao dos dfs comas supostas novas variaveis de interesse
imoveis_clean_rua <- imoveis_clean %>%
group_by(rua) %>%
summarise(total_a = n(),
media_p = mean(price_by_sqm, na.rm = T)) %>%
filter(!(is.na(rua)))
imoveis_clean_condominio <- imoveis_clean %>%
group_by(condominio) %>%
summarise(total_a = n(),
media_p = mean(price_by_sqm, na.rm = T)) %>%
filter(!(is.na(condominio)))
# como ha agencias sem nenhuma entrada, filtro antes de achar as variaves finais
imoveis_clean_agent <- imoveis_clean %>%
group_by(agent) %>%
summarise(total_a = n(),
media_p = mean(price_by_sqm, na.rm = T)) %>%
filter(!(is.na(agent))) %>%
filter(!(is.na(media_p)))
# plot total_a
var_t <- list(imoveis_clean_rua %>% arrange(desc(total_a)),
imoveis_clean_condominio %>% arrange(desc(total_a)),
imoveis_clean_agent %>% arrange(desc(total_a)))
plots_total <- list()
for (i in 1:length(var_t)){
plots_total[[i]] <- ggplot(var_t[[i]], aes(x = factor(var_t[[i]][[1]], levels = var_t[[i]][[1]]),
y = total_a)) +
geom_bar(stat="identity", width=.5) +
coord_flip() +
labs(title=paste('Anúncios x ', names(var_t[[i]])[[1]]),
x = names(names(var_t[[i]])[[1]]),
y = 'Anúncios')
print(plots_total[[i]])
}
# plot da media_p
var_m <- list(imoveis_clean_rua %>% arrange(desc(media_p)),
imoveis_clean_condominio %>% arrange(desc(media_p)),
imoveis_clean_agent %>% arrange(desc(media_p)))
plots_media <- list()
for (i in 1:length(var_m)){
plots_media[[i]] <- ggplot(var_m[[i]], aes(x = factor(var_m[[i]][[1]], levels = var_m[[i]][[1]]),
y = media_p)) +
geom_bar(stat="identity", width=.5) +
coord_flip() +
labs(title=paste('Media x ', names(var_m[[i]])[[1]]),
x = names(names(var_m[[i]])[[1]]),
y = 'Media')
print(plots_media[[i]])
}
# variaveis criadas
vars_rua <- imoveis_clean_rua %>% select(rua, anuncio_rua = total_a, valores_rua = media_p)
vars_condo <- imoveis_clean_condominio %>% select(condominio, anuncio_condo = total_a, valores_condo = media_p)
vars_agente <- imoveis_clean_agent %>% select(agent, anuncio_agente = total_a, valores_agente = media_p)
Após terminar a construção de novas variáveis, prossegue-se às etapas de criação do dataset final para análise. Começarei pela junção das novas variáveis, e exclusão das identificadoras, além de ‘bairro’ que, como visto ao longo de toda análise exploratória, possui variância quase inexistente.
# df final para analise e com variaveis novas
imoveis_final <- imoveis_clean %>%
left_join(roomsBathGarage) %>%
left_join(vars_rua) %>%
left_join(vars_condo) %>%
left_join(vars_agente) %>%
select(-c(Data, rua, numero, bairro, price_by_sqm, condominio, rooms, bathrooms, garage,agent, agent_number,
latitude, longitude, property_d, url))
## Joining, by = "property_d"
## Joining, by = "rua"
## Joining, by = "condominio"
## Joining, by = "agent"
glimpse(imoveis_final)
## Observations: 7,295
## Variables: 13
## $ price <int> 1000000, 1000000, 1000000, 1000000, 1000000, 1...
## $ area <int> 71, 67, 45, 79, 70, 70, 48, 85, 90, 100, 73, 7...
## $ condominium_fee <dbl> 1450, 1280, 700, 648, 1090, 1100, 1100, 800, 7...
## $ iptu <dbl> 280, 2200, NA, 220, 240, 240, 120, 140, 100, N...
## $ rooms_coerc <dbl> 2, 2, 1, 1, 1, 1, 1, 1, 1, 3, 2, 1, 2, 2, 3, 1...
## $ bathrooms_coerc <dbl> 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 2, 2, 2, 2, 4, 1...
## $ garage_coerc <dbl> 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1...
## $ anuncio_rua <int> 320, 320, 49, 570, 2, 94, 383, 276, 276, 383, ...
## $ valores_rua <dbl> 18.31858, 18.31858, 19.19827, 18.79077, 15.164...
## $ anuncio_condo <int> NA, 29, 7, 87, NA, 21, 85, 53, 53, 10, 18, NA,...
## $ valores_condo <dbl> NA, 15.17407, 22.18257, 19.48547, NA, 16.90986...
## $ anuncio_agente <int> 1, 3, 2, 72, 9, 73, 49, 49, 2, 19, 3, 5, 74, 1...
## $ valores_agente <dbl> 14.08400, 13.78067, 21.11100, 18.12724, 16.408...
Em seguida, prosseguir-se-á “fatorização” do dataframe, ou seja, transformação das variáveis categóricas para o formato factor
# Função mais eficiente para featurização
imoveis_fact <- imoveis_final %>%
mutate_at(
.vars = vars('rooms_coerc', 'bathrooms_coerc','garage_coerc'),
.funs = funs(as.factor(.))
)
glimpse(imoveis_fact)
## Observations: 7,295
## Variables: 13
## $ price <int> 1000000, 1000000, 1000000, 1000000, 1000000, 1...
## $ area <int> 71, 67, 45, 79, 70, 70, 48, 85, 90, 100, 73, 7...
## $ condominium_fee <dbl> 1450, 1280, 700, 648, 1090, 1100, 1100, 800, 7...
## $ iptu <dbl> 280, 2200, NA, 220, 240, 240, 120, 140, 100, N...
## $ rooms_coerc <fct> 2, 2, 1, 1, 1, 1, 1, 1, 1, 3, 2, 1, 2, 2, 3, 1...
## $ bathrooms_coerc <fct> 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 2, 2, 2, 2, 4, 1...
## $ garage_coerc <fct> 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1...
## $ anuncio_rua <int> 320, 320, 49, 570, 2, 94, 383, 276, 276, 383, ...
## $ valores_rua <dbl> 18.31858, 18.31858, 19.19827, 18.79077, 15.164...
## $ anuncio_condo <int> NA, 29, 7, 87, NA, 21, 85, 53, 53, 10, 18, NA,...
## $ valores_condo <dbl> NA, 15.17407, 22.18257, 19.48547, NA, 16.90986...
## $ anuncio_agente <int> 1, 3, 2, 72, 9, 73, 49, 49, 2, 19, 3, 5, 74, 1...
## $ valores_agente <dbl> 14.08400, 13.78067, 21.11100, 18.12724, 16.408...
Chega-se à parte mais relevante do trabalho: construção do modelo preditivo para precificação de imóveis. O pacote escolhido para a tarefa foi o “mlr”. Atualmente, um dos pacotes mais utilizados por usuários da linguagem, possui muitas opções de customização e facilidade de aplicação.
A primeira etapa consistirá na preparação dos dados através da padronização das variáveis contínuas, o tratamento dos valores missing e o processo conhecido como “One Hot Encoding”, ou seja, transformação das variáveis factor em dummies por levels.
## Pacote MLR
# Pacote
library(mlr)
## Loading required package: ParamHelpers
set.seed(666)
# Imputador dos valores missing
imp <- impute(
imoveis_fact,
classes = list(
factor = imputeMode(),
integer = imputeMean(),
numeric = imputeMean()
)
)
imoveis_imp <- imp$data
# Sumarizando as colunas
summarizeColumns(imoveis_imp) %>%
knitr::kable(digits = 2) # gerador de tabelas muito prático do knitr
| name | type | na | mean | disp | median | mad | min | max | nlevs |
|---|---|---|---|---|---|---|---|---|---|
| price | numeric | 0 | 4253932.96 | 4780029.51 | 2900000.00 | 2824353.00 | 250000.00 | 4.800e+07 | 0 |
| area | numeric | 0 | 211.11 | 158.71 | 168.00 | 133.43 | 18.00 | 1.683e+03 | 0 |
| condominium_fee | numeric | 0 | 2860.32 | 5790.93 | 2450.00 | 1868.08 | 1.00 | 4.600e+05 | 0 |
| iptu | numeric | 0 | 2073.65 | 3929.10 | 1500.00 | 1156.43 | 1.00 | 5.000e+04 | 0 |
| rooms_coerc | factor | 0 | NA | 0.63 | NA | NA | 1108.00 | 2.698e+03 | 4 |
| bathrooms_coerc | factor | 0 | NA | 0.72 | NA | NA | 899.00 | 2.034e+03 | 6 |
| garage_coerc | factor | 0 | NA | 0.73 | NA | NA | 1221.00 | 1.970e+03 | 5 |
| anuncio_rua | numeric | 0 | 245.57 | 137.24 | 245.57 | 113.52 | 1.00 | 5.700e+02 | 0 |
| valores_rua | numeric | 0 | 17.65 | 3.34 | 17.65 | 1.75 | 6.25 | 4.866e+01 | 0 |
| anuncio_condo | numeric | 0 | 38.53 | 22.09 | 38.53 | 0.00 | 1.00 | 1.180e+02 | 0 |
| valores_condo | numeric | 0 | 17.49 | 3.95 | 17.49 | 0.00 | 4.20 | 5.177e+01 | 0 |
| anuncio_agente | numeric | 0 | 42.37 | 33.57 | 38.00 | 37.06 | 1.00 | 1.350e+02 | 0 |
| valores_agente | numeric | 0 | 17.70 | 2.81 | 18.02 | 2.10 | 2.43 | 5.839e+01 | 0 |
# Normalização
imoveis_norm <- normalizeFeatures(imoveis_imp, target = "price")
# One hot encoding
imoveis_preProc <- createDummyFeatures(
imoveis_norm, target = "price",
cols = c(
"rooms_coerc",
"bathrooms_coerc",
"garage_coerc"
)
)
Com os dados preparados para análise, cria-se o objeto básico utilizado neste pacote: task (no caso, uma task específica para regressão). A partir dele, faz-se a separação hold-out (treino x teste)
# Criação do Task para Regressão (objeto usado nas operações do pacote)
task <- makeRegrTask(id = 'Imoveis_SP', data = imoveis_preProc, target = 'price')
# train test split (holdout é o método de resample de separar em train/test)
holdout <- makeResampleInstance('Holdout', task)
task_train <- subsetTask(task, holdout$train.inds[[1]])
task_test <- subsetTask(task, holdout$test.inds[[1]])
Com os objetos criados, pode-se iniciar o treinamento dos modelos. Serão usados 2 neste trabalho: um linear simples sem qualquer tipo de tratamento, para se usar como benchmark, e um XGBoost usando validações cruzadas e tunagem de hiper-parâmetros. A escolha do modelo deveu-se pelas relações não lineares entre as variáveis númericas, e seu sistema de uso de gradiente para minimização dos resíduos. 5 validações cruzadas serão utilizadas para evitar o overfitting e a série de hiper-parâmetros escolhidos para tunagem podem ter sua explicação melhor desenvolvida em https://github.com/dmlc/xgboost.
## Modelo linear básico
# modelo padrão linear
regr.lm <- makeLearner(id = 'lm', 'regr.lm')
# treinamento modelo
lm <- train(regr.lm, task_train)
## XGBoost
# Criação do modelo simples
xgb_model <- makeLearner("regr.xgboost")
## Warning in makeParam(id = id, type = "numeric", learner.param = TRUE, lower = lower, : NA used as a default value for learner parameter missing.
## ParamHelpers uses NA as a special value for dependent parameters.
# Grid de alguns parâmetros
xgb_params <- makeParamSet(
# The number of trees in the model (each one built sequentially)
makeIntegerParam("nrounds", lower = 100, upper = 2000),
# number of splits in each tree
makeIntegerParam("max_depth", lower = 1, upper = 10),
# "shrinkage" - prevents overfitting
makeNumericParam("eta", lower = .1, upper = .5),
# L2 regularization - prevents overfitting
makeNumericParam("lambda", lower = -1, upper = 0, trafo = function(x) 10^x)
)
# Controle da tunagem randômica
control <- makeTuneControlRandom(maxit = 1) # maxit represente tempo máximo de 1 min para iterações
# Plano de resampling
resample_desc <- makeResampleDesc("CV", iters = 5)
# tunador dos hiper parâmetros
tuned_params <- tuneParams(
learner = xgb_model,
task = task_train,
resampling = resample_desc,
measures = rmse, # R-Squared performance measure, this can be changed to one or many
par.set = xgb_params,
control = control
)
## [Tune] Started tuning learner regr.xgboost for parameter set:
## Type len Def Constr Req Tunable Trafo
## nrounds integer - - 100 to 2e+03 - TRUE -
## max_depth integer - - 1 to 10 - TRUE -
## eta numeric - - 0.1 to 0.5 - TRUE -
## lambda numeric - - -1 to 0 - TRUE Y
## With control class: TuneControlRandom
## Imputation value: Inf
## [Tune-x] 1: nrounds=1774; max_depth=3; eta=0.465; lambda=0.374
## [Tune-y] 1: rmse.test.rmse=1209875.5930818; time: 0.3 min
## [Tune] Result: nrounds=1774; max_depth=3; eta=0.465; lambda=0.374 : rmse.test.rmse=1209875.5930818
# Modelo tunado
xgb_tuned_model <- setHyperPars(
learner = xgb_model,
par.vals = tuned_params$x
)
# Verificação da performance do modelo tunado usando as validações cruzadas usadas para tunagem
resample(xgb_tuned_model,task_train,resample_desc,measures = list(rmse, rsq))
## Resampling: cross-validation
## Measures: rmse rsq
## [Resample] iter 1: 1009558.34928650.9477924
## [Resample] iter 2: 1616862.04355040.8780969
## [Resample] iter 3: 1167638.40041130.9387638
## [Resample] iter 4: 1076768.89157160.9451954
## [Resample] iter 5: 1423415.55945190.9258245
##
## Aggregated Result: rmse.test.rmse=1279247.7057512,rsq.test.mean=0.9271346
##
## Resample Result
## Task: Imoveis_SP
## Learner: regr.xgboost
## Aggr perf: rmse.test.rmse=1279247.7057512,rsq.test.mean=0.9271346
## Runtime: 20.0695
# Treinamento
XGBoost <- train(xgb_tuned_model, task_train)
Agora, com os modelos devidamente treinadss, podemos fazer as comparações dos resultados, quais variáveis possuem maior efeito e como é a sua influência nas previsões. A começar pela performance nos dados de teste, faz-se uma comparação de 2 indicadores de qualidade das previsões juntamente a uma ilustrativa.
# predição nos dados de teste
pred_lm <- predict(lm, task_test)
## Warning in predict.lm(.model$learner.model, newdata = .newdata, se.fit =
## FALSE, : prediction from a rank-deficient fit may be misleading
performance(pred_lm, measures = list(rmse, rsq))
## rmse rsq
## 3.124175e+06 5.905637e-01
pred_xgb <- predict(XGBoost, task_test)
performance(pred_xgb, measures = list(rmse, rsq))
## rmse rsq
## 1.196479e+06 9.399483e-01
## Análise dos dados fitted e feature importance
comparacao <- data_frame(pred_lm = pred_lm[['data']][['response']], pred_xgb = pred_xgb[['data']][['response']],
price = getTaskData(task_test)[['price']])
# Plot predictions (on x axis) vs actual bike rental count
comparacao %>%
tidyr::gather(tipo, valor, price, pred_lm, pred_xgb) %>%
filter(!(tipo == 'pred_xgb')) %>%
ggplot(aes(valor, fill = tipo)) +
ggtitle('LM') +
geom_density(alpha = .5)
comparacao %>%
tidyr::gather(tipo, valor, price, pred_lm, pred_xgb) %>%
filter(!(tipo == 'pred_lm')) %>%
ggplot(aes(valor, fill = tipo)) +
ggtitle('LM') +
geom_density(alpha = .5)
Como demonstrado acima, e já esperado, o modelo xgbost apresentou índices superiores nos indicadores, além de, visualmente, ter uma aderência maior aos dados de teste. Portnato, para encerramento da análise, usarei o pacote ‘iml’, recomendado pelos criadores do ‘mlr’, para aumentar a compreensão desse modelo mais complexo. Primeiro, descobrindo qual a importância de cada uma das features e, em seguida, como os valores da mais importante influenciam os valores previstos pelo modelo.
# Pacote
library(iml)
# usando Predictor$new() para criar um plot das features mais importantes
X = getTaskData(task_train)[getTaskFeatureNames(task_train)]
predictor = Predictor$new(XGBoost, data = X, y = getTaskData(task_train)['price'])
imp = FeatureImp$new(predictor, loss = "rmse")
plot(imp)
# Plot de como a variável mais importante afeta as previsões na média
pdp_area = Partial$new(predictor, feature = "area")
## Warning: The FeatureEffect class replaces the Partial class. Partial will
## be removed in future versions.
plot(pdp_area)